人类胸腺发育的细胞图谱揭示了T细胞组库的形成

分享是一种态度

文章信息

文献标题：A cell atlas of human thymic development defines T cell repertoire formation
发表时间：2020.02.21
发表杂志：Science（IF=41.845）
原文链接：http://dx.doi.org/10.1126/science.aay3224

摘要

INTRODUCTION：胸腺（thymus）是 T 细胞发育和 T 细胞受体（T cell receptor, TCR）组库形成的重要器官，塑造了机体的适应性免疫。T 细胞的胸腺内发育有空间协调性，受胸腺微环境（thymic microenvironment）多种细胞类型的精细调节。尽管胸腺在多种动物模型中被广泛研究，目前尚缺少一份完整的人类胸腺图谱帮助我们理解人体免疫系统。

RATIONALE：作者将处于发育中、儿童期及成人期的胸腺标本解离后，利用单细胞转录组测序构建了跨越人类生命周期的完整胸腺图谱。作者总共采集了 15 个胎儿胸腺（从 7 - 17 孕周）和 9 个出生后至成人期的胸腺样本，同时采用了不同的细胞分选策略保证细胞亚群的覆盖率。根据 scRNA-seq 得到的 marker 基因，作者利用单分子荧光原位杂交（single-molecule fluorescence in situ hybridization, smFISH）对不同状态的细胞作了空间定位。另外，作者系统地比较了人和小鼠胸腺单细胞转录组数据的异同。最后，作者还获取了单细胞 TCR 组库数据以探索人类 TCR 重排和选择的时空偏倚。

RESULTS：作者从人类胸腺中鉴定出了超过 50 种状态的细胞。人胸腺中的细胞状态、丰度及基因表达谱随着胎儿及出生后发育呈动态变化。作者鉴定出了新的胸腺成纤维细胞和上皮细胞，并对其进行了定位。作者通过计算预测了人 T 细胞从胎儿肝脏造血干祖细胞到多种成熟 T 细胞类型的发育轨迹，进而构建了导致 T 细胞命运决定的转录因子网络框架。在非经典的胸腺 T 细胞中，作者注意到了一种 CD8αα⁺ T，以表达 GNG4 为特征，定位于胸腺髓质周围，它们高表达 XCL1，并且与 XCR⁺ 树突状细胞共定位。与小鼠数据相比，作者发现人体中这类非经典胸腺 T 细胞有着不同的基因表达谱。最后，作者发现人 VDJ 在重排和多轮选择下呈现出强烈偏倚，包括 CD8⁺ T 细胞中的 TCRα V-J 偏倚等。

CONCLUSION：该文章描绘了跨生命周期、跨物种的胸腺单细胞转录组图谱，提供了胸腺微环境下 T 细胞发育的高分辨率景观。通过系统地跨物种比较，作者揭示了人类特异性的细胞状态和基因表达特征。这份详细的细胞网络将有助于未来构建体外类器官模型，用以还原人类的在体胸腺组织。

测序数据

• 组织样本：

1. 人胸腺样本：胚胎（7 - 17 PCW）、出生后（3 m - 15 y）、成人（25 y、35 y）
2. 小鼠胸腺样本：C57BL/6J （4、8、24 w）

分选策略：

1. DAPI 阴选富集活细胞
2. CD45 分别富集免疫细胞和非免疫细胞
3. CD3 分别富集 T 细胞和非 T 细胞
4. EPCAM 富集上皮细胞

建库方法：

1. 10x Genomics Single Cell 3' v2 and 5’ Reagent Kit：人类样本
2. 10x Genomics Single Cell 3' v3 Reagent Kit：小鼠出生后样本
3. 10x Genomics Single Cell 5' VDJ：部分人类样本
4. Smart-seq2：部分人类样本

公共数据：

1. scRNA-seq：人胚胎肝脏（E-MTAB-7407）、人骨髓造血干细胞（GSE117498）、小鼠胚胎胸腺（GSE107910）、小鼠胸腺基质细胞（GSE103967）
2. Microarray：人出生后分选的胸腺细胞（E-MEXP-337）

数据分析

由于本文大量篇幅用在阐述 10x 数据，Smart-seq2 数据主要用于验证，这里我重点总结 10x 数据的分析方法。

原始数据处理

基因表达数据：Cell Ranger Single-Cell Software Suite (version 2.0.2 for 3’ v2 chemistry, version 3.0.0 for 3’ v3 chemistry and version 2.1.0 for 5’ chemistry)

VDJ 数据：Cell Ranger Single-Cell Software Suite (version 2.1.0)

下游分析

基本分析流程：

python3 下使用 scanpy 软件（version 1.3.4）

1. 归一化：scanpy.api.pp.normalize_per_cell
2. log 转换：scanpy.api.pp.log1p
3. 高可变基因（Highly variable genes, HVGs）：scanpy.api.pp.filter_gene_dispersion
4. 细胞周期：将预先定义的细胞周期相关基因从 HVGs 中移除
5. Scaling data：scanpy.api.pp.scale
6. PCA：scanpy.api.pp.pca 基于高可变基因
7. 数据整合：scanpy.api.pp.bbknn
8. 聚类：scanpy.api.tl.louvain

质控：

1. 保留 UMI > 2000、500 < Gene < 7000 的细胞
2. 鉴定 doublets：scrublet（参数：sim_doublet_ratio=2、n_neighbors=30、 expected_doublet_rate= 0.1），scrublet 评分 > 0.7 为 doublet
3. 过度预聚类（sc.tl.louvain 分辨率 20），平均 scrublet 评分 > 0.6 的聚类为 doublets。剩余的聚类中，若符合以下标准则同样推断为 doublets：
   (1) 根据先验知识，同时表达两种或以上细胞类型 marker 基因（例如同时表达 CD3 和 CD19）
   (2) UMI 数过高
   (3) 缺少定义该聚类的 marker 基因
4. 移除甲状腺来源的污染细胞聚类：Epi_PAX8（PAX8, HHEX, TG, NKX2.1）和 Epi_GCM2（PTH, GCM2, GATA3, CHGA）

批次效应矫正：

考虑到批次效应可能来自技术差异（例如 10x 试剂版本）和生物学差异（例如发育阶段和组织供体），作者采用了迭代整合方式：

1. 利用 scanpy.api.pp.bbknn 整合相似样本（例如所有胎儿和儿童来源样本）
2. 利用整合后的数据获取粗略的细胞亚群注释
3. 利用批次信息（如 10x 试剂版本、供体等）或细胞亚群注释作为分类变量，拟合 L2-正则化的线性模型（L2-regularised linear model），将批次信息造成的差异回归掉（regress out），保留仅包含生物学信息的残差
4. 再次利用 scanpy.api.pp.bbknn 整合数据，用于下游的细化聚类、可视化和轨迹推断

聚类和细胞亚群注释：

1. 根据已知 marker 基因和差异表达基因（differentially expressed genes, DEGs）进行手动注释，首先将分群清晰、具有唯一定义的细胞聚类注释出来
2. 利用上述注释的聚类训练 logistic 回归模型（logistic regression model），用以预测因转录相似性而被软件聚到一起的混合细胞亚群
3. 对粗略注释的细胞分群（如上皮细胞、单阳性 T 细胞等），将其提取出来重复 HVGs 选择、降维和聚类流程，以获取更精细的细胞亚群（如髓质胸腺上皮细胞 mTEC 或调节性 T 细胞等）

估计细胞类型的比例：

1. 作者首先宽泛地定义细胞类型（如淋巴细胞、髓细胞等），并计算每种类型在选定的对比组之间的比例
2. 如果一次比较中的所有细胞类型均来自同一个分选门，则将比例简单定义为：特定类型细胞数 / 细胞总数
3. 如果细胞来自不同的分选门，则为每个分选门计算一个归一化因子：给定分选门的细胞数 / 所有分选门的总细胞数。对每个分选门的细胞数均乘以相应的归一化因子，再用归一化的细胞数来计算细胞比例
4. 用 t 检验评估细胞比例变化的显著性

轨迹推断：

作者利用前述批次矫正方法获得邻域图（neighbourhood graph），选择符合先验知识（包括已知 T 细胞分化、TCR 重排的顺序）的流形（manifold），然后利用 scanpy.api.tl.dpt 计算扩散拟时间（diffusion pseudotime），并计算拟时间中差异表达的基因

转录因子网络：

1. 人类转录因子信息下载自数据库 AnimalTFDB3
2. 基因表达估算：在 3D UMAP 空间对最近 30 邻域（30-nearest neighbors）取平均
3. 在 3D UMAP 的每个体素中随机抽样细胞，以去除冗余信息。在用于估算表达量的最近 30 邻域内计算细胞类型频率，对每个细胞类型计算注释评分
4. 选取高可变的转录因子计算相关性矩阵，建立图空间，并利用力导引图（force-directed graph, FDG）实现可视化

• VDJ 序列分析：

1. 质控：选取包含全长的重排序列且 UMI 数 > 2 的 TCR 链
2. 鉴定生产性重排（productive）的 TCR：能找到覆盖互补决定区（complementarity determining regions, CDR）CDR3 的开放读码框（Open Reading Frame, ORF）
3. 比较 V、D、J 基因的使用率：对每个细胞类型的每个 VDJ 基因计数，利用该细胞类型的总 VDJ 计数进行归一化，并转换为 z-score，利用 t 检验在不同细胞类型之间比较 z-score。对 CD4⁺ T 和 CD8⁺ T 细胞的比较，还使用了 Cochran–Mantel–Haenszel 检验

细胞间相互作用：CellPhoneDB(www.CellPhoneDB.org)

1. 每种细胞类型抽样 1000 个细胞
2. 从 HUGO Gene Nomenclature Committee 数据库获取趋化因子信息
3. 通过在细胞-细胞对内将配-受体对的平均表达水平相乘，计算相互作用评分，并使用最大值归一化
4. 根据特异性的显著程度选取相互作用对（interaction pairs）进行可视化

整合公共数据：

同样使用了 BBKNN 方法

跨物种比较：

1. 使用前述流程对小鼠数据集做预处理，包括归一化、批次整合、降维聚类等。利用人类数据集训练 logistic 模型对小鼠数据进行辅助注释
2. 数据整合：结合了两种方式：
   (1) 基于互为训练集的 logistic 模型做双向预测，两套预测概率相乘作为最终的相似度评分
   (2) 线性回归（linear regression）辅助的 BBKNN 方法去除物种特异性变异，保留生物学信息的结构

主要结果

人类生命过程中胸腺内的细胞组成

• 质控后保留了 255,901 个单细胞，包含了从胎儿期、幼年期、青春期到成人期的完整胸腺图谱
• 鉴定出超过 40 种细胞亚型：
• 免疫细胞：包括分化的 T 细胞（DN、DP、CD4 / CD8 SP、T_reg 、CD8αα⁺ T 和 γδ T）、B 细胞、NK 细胞、固有免疫细胞（ILCs）、巨噬细胞、单核细胞和树突状细胞（dendritic cells, DCs）等
• 非免疫细胞：胸腺上皮细胞（thymic epithelial cells, TECs）、成纤维细胞、血管平滑肌细胞（vascular smooth muscle cells, VSMCs）、内皮细胞和淋巴管内皮细胞等
• DCs 进一步被分为经典的髓样 DCs（DC1、DC2）和浆细胞样 DC（plasmacytoid DC, pDC）
• 成纤维细胞细被细分为 Fb1（COLEC11, C7, GDF10）、Fb2（PI16, FN1, FBN1）和  Fb_cycling。Fb1 表达固有免疫相关的重要基因 COLEC11 以及调节上皮细胞发育的维甲酸反应酶 ALDH1A2。而 Fb2 表达细胞外基质基因和信号素（semaphorins），调控血管发育。smFISH 发现 Fb1 主要定位在胸腺小叶周围，而 Fb2 主要在小叶间区，并且和大血管相联系。
• 作者将 TECs 根据髓质（medullary）和皮质（cortical）分为 mTEC 和 cTEC，并利用小鼠数据集辅助注释，鉴定出 cTEC（PSMB11）、mTEC(I)（KRT14）、mTEC(II)（AIRE）、mTEC(III)（KRT1）以及罕见的 mTEC(IV)（DCLK1, POU2F3）。作者发现 cTECs 在早期发育阶段丰度较高，而在胎儿晚期和儿童期则出现一种中间态亚型 mcTEC（DLK2）。
• 作者还注意到两群 EpCAM⁺ 细胞：肌样 TEC（MYOD1, MYOG）和神经内分泌样 TEC（NEUROD1, NEUROG1, CHGA），定位于胸腺髓质。自身免疫性重症肌无力相关的基因 CHRNA1 在这两群细胞以及 mTEC(II) 中高表达。
• 最后，作者分析了已知的先天性 T 细胞免疫缺陷（congenital T cell immunodeficiencies）致病基因表达谱，揭示了这些基因可能发挥作用的位置和发育阶段

胸腺基质和 T 细胞的协调发育

• 早期胎儿胸腺（7 - 8 PCW）中，淋巴细胞主要包括 NK 细胞、γδ T 细胞和 ILC3；分化中的 αβ T 很少，主要处于 DN 阶段。随后 DN 逐渐分化出 DP 和 SP 阶段，并在 12 PCW 达到均衡；相反，固有淋巴细胞比例则逐渐减少
• 逐渐退化的成人胸腺中出现终末分化的 T 细胞，可能是来自归巢 T 细胞或外周血 T 细胞。退化胸腺中，细胞毒性 CD4⁺ T 细胞（CD4⁺ CTLs）、记忆 T 细胞和记忆 B 细胞比例均增加
• T 细胞发育的趋势和胸腺基质细胞相对应。随着 T 细胞成熟的启动，TEC 从 cTEC 主导逐渐向 cTEC 和 mTEC 平衡状态转变，提示了胸腺上皮细胞和成熟 T 细胞相互作用、同步分化的特征
• 成纤维细胞中，Fb1 在早期发育阶段占主导，而后 Fb1 和 Fb2 逐渐达到均衡
• 其他免疫细胞中，巨噬细胞在早期发育中较丰富，而 DCs 随着发育进程逐渐增加。其中，DC1 在 12 PCW 后占主导，而 pDC 在出生后迅速增加
• 作者最后预测了胸腺基质细胞和 T 细胞相互作用的配-受体对：淋巴毒性信号（LTB:LTBR）来自多种免疫细胞，并由大部分基质细胞接受；而 RANKL-RANK 信号（TNFRSF11:TNFRSF11A）局限于 ILC3 和 mTEC(II) / 淋巴管内皮细胞；FGF 信号（FGF7:FGFR2）从成纤维细胞到 TEC，且 FGFR2 在成人胸腺表达降低；NOTCH1 是早期胸腺祖细胞（early thymic progenitors, ETPs）的主要受体，而多种 Notch 配体中，JAG2 和 DLL4 主要由 cTECs 和内皮细胞表达，其他 TEC 则广泛表达 JAG1

传统 T 细胞分化轨迹

• 作者整合了胎儿肝脏造血干祖细胞和胸腺 T 细胞的数据，发现 ETP 在 UMAP 图上位于造血干细胞 / 多能干祖细胞（hematopoietic stem cells / multipotent progenitors, HSCs/MPPs）和 pre / pro-B 细胞之间。而进一步整合骨髓单细胞数据后，作者发现 ETP 位于骨髓多重淋巴祖细胞（multi-lymphoid progenitor, MLP）和胎儿肝脏早期淋巴祖细胞（early lymphoid progenitor）之间

• 作者进一步推断了 T 细胞后续分化的轨迹：始于 CD4^-CD8^- DN，逐渐变为 CD4⁺CD8⁺ DP，经过 CCR9 ^high 的 αβ(entry) 阶段，分化为成熟的 CD4⁺ 或 CD8⁺ SP。同时，作者注意到 γδ T 在 DN-DP 交界处另行分化。

• 对 DN 和 DP，作者根据其细胞周期增殖（CDK1）将其细分为增殖态（proliferating）和静息态（quiescent）。VDJ 重排基因（RAG1, RAG2）从增殖晚期开始增加，于静息态达到高峰，提示了 T 细胞分化中增殖和重排的关系

• 通过整合 TCR 数据，作者发现，在增殖晚期 DN 中检测到重排的 TCRβ 链，与重排事件和 pre-TCRα（PTCRA）表达的增加相符合。TCRβ 链非生产性 / 生产性重排的比率在 DN 中相对较高，而在进入 DP 阶段逐渐下降到基线水平，提示了 β-selection 的影响。在 DN(Q) 阶段，TCRβ 链非生产性重排达到高峰，可能提示第一轮等位基因重排失败的细胞进行了其他等位基因重排。在 DP 阶段，自增殖期开始检测到重排的 TCRα 链，且与 TCRβ 链不同，其在 DP(Q) 阶段非但没有富集，反而直接被丢弃了

• 通过推断拟时间差异表达基因，作者鉴定了新的 T 细胞发育阶段及 marker 基因，例如early DN（ST18）、DP （AQP3）和 DP-to-SP 过渡态（TOX2），并且构建了 T 细胞谱系定向相关的转录因子调控网络

• 最后，作者结合基于蛋白 marker 分选的 microarray 公共数据，发现 DN(P)、DN(Q) 和 DP(P) 分别对应 CD34⁺CD1A⁺ DN、CD4⁺ ISP 和 CD3^- DP，而 DP(Q) 和 αβ(entry) 对应 CD3⁺ DP。其中 DN(Q) 和 CD4⁺ ISP 均表达 pre-β-selection 特征

T_reg 和非经典 T 细胞的发育

• 作者鉴定出了多种非经典 T 细胞，例如 T_reg（FOXP3）、CD8αα⁺ T、NKT 样细胞（EOMES, NKG7, IFNG, TBX21）以及 T_H17 样细胞（CD4, CD40LG, RORC,  CCR6）。这些细胞的发育需要激动剂选择（agonist selection），其非生产性 TCR 链比例很低，提示在胸腺内的时间比经典 T 细胞更久
• 非经典 T 细胞在胸腺的富集程度高于肝脏，并且均富集于成熟后（~10 PCW）的胸腺，提示其发育具有胸腺依赖性
• 在 UMAP 图中，T_reg 和 αβ T 有分化轨迹联系。作者将连接处的 T_reg 定义为 T_reg(diff)，其表达低水平 FOXP3 和 CTLA4，而高表达 IKZF4、GNG8 和 PTGIR，这些基因均与自身免疫和 T_reg 分化相关。此外，作者注意到另一群由非编码 RNA（MIR155HG）定义的不同于终末分化 T_reg 及 T_reg(diff) 的细胞，注释为 T_(agonist) ，其表达 IL2RA 但低表达 FOXP3，与小鼠胸腺中 CD25⁺FOXP3^- T_reg 祖细胞相似。T_reg(diff) 和 T_(agonist) 可能代表了人类胸腺中的两种 T_reg 祖细胞
• 作者进一步将 CD8αα⁺ T 细胞分为 3 群：CD8αα⁺ T(I)（GNG4, CREB3L3, CD72）、CD8αα⁺ T(II)（ZNF683, MME）和 CD8αα⁺ NKT 样细胞（EOMES）。CD8αα⁺ T(I) 和 CD8αα⁺ T(II) 在早期高表达 PDCD1，而在终末分化阶段逐渐下调。CD8αα⁺ T(I) 和晚期 DP 及 αβ(entry) 之间有明显的谱系区分，CD8αα⁺ T(II) 则混合了 αβ 和 γδ T 的特征，而 NKT 样细胞则富集 γδ TCR 链
• 利用 smFISH，作者发现 GNG4⁺ CD8αα⁺ T(I) 细胞位于胸腺髓质。CD8αα⁺ T(I) 和 T_reg 都表达表面蛋白 marker CD137（TNFRSF9），于是作者通过流式分选和 Smart-seq2 建库测序，验证了 CD3⁺CD137⁺CD4^- 可以作为 CD8αα⁺ T(I) 的分选策略。跨物种比较发现，人 CD8αα⁺ T(I) 与小鼠上皮内淋巴细胞前体 A 型（intraepithelial lymphocyte precursor type A, IELpA）最为相似，均表达 HIVEP3、NR4A3、PDCD1 和 TNFRSF9，但同时也存在物种间差异基因（人类 GNG4 和 XCL1，小鼠 ZEB2 和 CLDN10），并且人类 CD8αα+ T(I) 最终成为 CD8A^highCD8B^low 表型，而小鼠 IELpA 则表现为 CD8A^lowCD8B^lowCD4^low 特征
• NKT 样细胞和 T_H17 样细胞均表达 KLRB1 和 ZBTB16，表现出类似固有淋巴细胞的特征
• 最后，作者分析了 TRDV1 和 TRDV2 两个常见的 TCRδ V 基因，发现 TRDV2 主要用于 DN 阶段，而 TRDV2 主要用于 DP(Q) 和 αβT entry 阶段。进一步地，作者推断 CD8αα+ T(II) 主要来自晚期 DP，而 NKT 样和 T_H17 样细胞主要来自早期分化的胸腺细胞

DCs 在胸腺细胞选择中的募集和活化

• DC 亚群及 marker 基因：前文所述的 DC1（XCR1, CLEC9A）、DC2（SIRPA, CLEC10A）和 pDC (IL3RA, CLEC4C），以及新鉴定的活化 DC（activated DCs, aDCs; LAMP3, CCR7）。aDC 高表达多种趋化因子、共刺激分子和转录因子（AIRE, FOXD4），它们可能代表了过去在人胸腺和扁桃体描述过的 AIRE⁺CCR7⁺ DCs
• aDC 进一步细分为 aDC1、aDC2 和 aDC3。作者发现 aDC1 与 aDC2 分别同 DC1 及 DC2 有相似表达谱，提示它们之间的激活关系（aDC1 - DC1 和 aDC2 - DC2），并且两群 DC 分别表达不同的趋化因子。aDC3 细胞中，主要组织相容性复合物（major histocompatibility complex, MHC）和共刺激分子的表达均相对降低，提示其可能为 DC 的活化后状态
• 作者重点关注了趋化因子介导的细胞间相互作用，发现 aDC 同时表达 CCR7 和 CCL19，提示其具有募集 T 细胞进入胸腺髓质的能力。aDC 还表达 CCL17 和 CCL22，相应的受体 CCR4 则由 CD4⁺ T 和某些特定的 T_reg 表达。aDC 还可能通过 CXCL9/10:CXCR3 募集其他 DC 和成熟 T_reg，并在 T_reg  成熟过程中提供共刺激信号。作者还注意到 CD8αα⁺ T(I) 表达 XCL1，可能和 DC1（XCR1）的募集与活化有关
• 最后，作者利用 smFISH 验证了上述发现

人 TCR 组库形成和选择的偏倚

• TCRβ 的 VDJ 基因重排偏倚从 DN 到成熟 T 细胞的过程中持续存在，该偏倚可能和基因座的襻环结构有关。D2 基因与 J2 基因之间具有偏好性，而 D1 基因和 J1、J2 重排的频率差异不大。V - D 或 V - J 之间则没有明显相关性。在 β-selection 后，某些 TRBV 基因被丢弃或富集，提示在不同 Vβ 基因对 MHC 肽刺激的反应能力具有种系编码的差异
• TCRα 基因座中，发育时间和 V - J 配对之间有清晰的相关性：近端先于远端重排，V - J 基因严格配对。而在成熟 T 细胞中，远端 V - J 配对被丢弃，可能反映了阳性选择（positive selection）阶段的偏倚
• 作者比较了不同细胞亚群之间的 TCR 组库差异，发现 CD8⁺ T 细胞和其他细胞有明显差别，其 TRAV - TRAJ 组库偏好远端 V - J 对，提示这种差异可能与 CD8⁺ T 谱系的缓慢或低效定向有关。和胚胎期相比，出生后胸腺中的 CD8αα⁺ T(I) 则偏好近端的 V - J 对

总结

本文构建了人类经典和非经典 T 细胞的分化轨迹以及 TCR 组库信息，发现了 TCR 组库在成熟 T 细胞中的偏倚，可能提示了人体如何应对不同的抗原挑战。本文对胸腺微环境的分析揭示了构成胸腺的复杂细胞类型，以及基质细胞和固有免疫细胞如何相互作用以支持 T 细胞分化。这些细胞间交流网络可以帮助体外构建产生 T 细胞的培养系统，并影响未来 T 细胞治疗的工程学策略。

我的评价

本文的生物学意义已经说得很清楚了，我来从生信分析的角度谈一些感想吧。

面对这样大量的单细胞样本，无论是数据读写还是内存占用，光靠基于 R 语言的分析平台已经远远不够。人类发育细胞图谱（HDCA）项目里，过去两篇母-胎界面（70,000 细胞）和胎儿肝脏造血（大约 130,000 细胞）的 Nature 基本还是以 Seurat 包为主，本文可能是该项目第一篇纯 python 分析的大文章吧。本文的代码已经以 jupyter notebook 的形式上传 github，包含大量经典套路和高级分析，其中 logistic 模型辅助细胞注释的部分感觉是 Teichmann 团队后来发表的 SCCAF(https://github.com/SCCAF) 方法的前身。因此，本文理论上说可以作为学习 python 环境下分析单细胞数据的案例。

于是乎，我们又来到了“可复现的生物信息学分析”这个话题。前两年，一位中科院计算所的老师曾告诉我：”即使公开了代码，生信分析的文章仍有 70 - 80 % 都不能完全重复出来。“然而，早在 2017 年，《Reproducible bioinformatics project: a community for reproducible bioinformatics analysis pipelines》(https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-018-2296-x)一文已经指出以 docker 为代表的容器化技术可能为我们带来更灵活可控的生信分析环境。理论上说，使用给定镜像的容器在任何设备上做分析，只要喂给计算机相同的数据，理应得出相同的结果。除此之外，各大杂志也在利用 github、codeocean 等代码托管或云计算平台推进生信分析规范化的进程，上文提到的 SCCAF作者同样提供了 docker 镜像用于重复该软件的结果。综上，未来的生信分析文章一定会面临更加严苛的检验。

话说我真的不是给 docker 打广告啊。醒醒，copy一下教程、随便跑跑代码出一堆图、水得一手文章的日子快要到头啦！

---

---

如果你对单细胞转录组研究感兴趣，但又不知道如何入门，也许你可以关注一下下面的课程

看完记得顺手点个“在看”哦！

生物 | 单细胞 | 转录组丨资料每天都精彩

长按扫码可关注